Resonant and Non-Resonant Modulated Amplitude Waves for Binary Bose-Einstein 

Condensates in Optical Lattices 



O 



> 



X 



Mason A. Porter^, P.G. Kevrekidis^, and Boris A. Malomed'^ 
^School of Mathematics and Center for Nonlinear Science, 
Georgia Institute of Technology, Atlanta GA 30332, USA 
^Department of Mathematics and Statistics, 
University of Massachusetts, Amherst MA 01003-4515, USA 
^Department of Interdisciplinary Studies, Faculty of Engineering, 
Tel Aviv University, Tel Aviv 69978, Israel 
(February 8, 2008) 



We consider a system of two Gross-Pitaevskii (GP) equations, in the presence of an optical-lattice 
' (OL) potential, coupled by both nonlinear and linear terms. This system describes a Bose-Einstein 

I condensate (BEG) composed of two different spin states of the same atomic species, which interact 

linearly through a resonant electromagnetic field. In the absence of the OL, we find plane-wave 
solutions and examine their stability. In the presence of the OL, we derive a system of ampli- 
tude equations for spatially modulated states which are coupled to the periodic potential through 
the lowest-order subharmonic resonance. We determine this averaged system's equilibria, which 
represent spatially periodic solutions, and subsequently examine the stability of the corresponding 
solutions with direct simulations of the coupled GP equations. We find that symmetric (equal- 
amplitude) and asymmetric (unequal-amplitude) dual-mode resonant states are, respectively, stable 
I and unstable. The unstable states generate periodic oscillations between the two condensate com- 

ponents, which is possible only because of the linear coupling between them. We also find four-mode 
states, but they are always unstable. Finally, we briefly consider ternary (three-component) con- 
Ch ', densates. 

PACS: 05.45.-a, 03.75.Lm, 05.30.Jp, 05.45.Ac 



I. INTRODUCTION 



> 
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' At sufficiently low temperatures, particles in a dilute boson gas condense in the ground state, forming a Bose- 
, Einstein condensate (BEG). This was first observed experimentally in 1995 in Na and Rb vapors [37,16,26,12]. 

In the mean-field approximation, a dilute BEG is described by the nonlinear Schrodinger (NLS) equation with an 
, external potential, which is also called the Gross-Pitaevskii (GP) equation. In particular, BEGs may be considered 
' in the quasi-one-dimensional (quasi-lD) regime, with the transverse dimensions of the condensate on the order of its 
mean healing length x (given by = (87rn|a|)~^) and a much larger longitudinal dimension [9,10,8,16]. The length 
i-H ' X is determined by the mean atomic density n and the two-body s-wave scattering length a, where the interactions 
^ , between atoms are repulsive if a > and attractive if a < [37,16,27,4]. 

The quasi-lD regime, which corresponds to "cigar-shaped" BEGs, is described by the ID limit of the 3D mean-field 
theory (rather than by a ID mean-field theory proper, which would only be appropriate for extremely small transverse 
dimensions of order ~ a) [9,8,10,45,6]. In this situation, the condensate wave function tp{x,t) obeys the effective ID 



GP equation, 



^2 

ihipt = -—ipxx + glipfi' + V{x)^p, 
2m 



where m is the atomic mass, V{x) is an external potential, g ~ [AnU a/m][l + 0(77^)], and rj = y/n\a\^ is 
the dilute-gas parameter [16,27,4,29]. Experimentally relevant potentials V{x) include harmonic traps and peri- 
odic potentials (created as optical lattices, which are denoted OLs and arise as interference patterns produced 
by coherent counterpropagating laser beams illuminating the condensate). In the presence of both potentials, 
V{x) = Vb cos [2^(2; — xo)] + Vix^/2, where xq is the offset of the periodic potential relative to the center of the 
the parabolic trap. When (2tt/k)'^ Vi <^Vo, the potential is dominated by its periodic component over many periods 
[17,14,11]; for example, when Vq/Vi — 500 and k — 10, the parabolic component in V{x) is negligible for the 10 
periods closest to the trap's center. In this work, we set Vi = and focus entirely on OL potentials. This assumption 
is motivated by numerous recent experimental studies of BEGs in OLs [22,3] and is widely adopted in theoretical 
studies [9,8,10,7,15,14,17,34,2,48,38,39,35,50,31,30]. 
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Multiple-component BECs, which constitute the subject of this work, have been considered in a number of the- 
oretical works [24,41,40,20,13,47,18]. Mixtures of two different species (such as *^Rb and *''Rb) are described by 
nonlinearly coupled GP equations: 




(1) 



(2) 



depends on the cross-scattering length ai2 [41]. There are numerous subcases of Eqs. (1) to consider, as various 
combinations of signs for the scattering coefficients gi, g2, and h may occur. It is important to note, however, that 
if 9i,2 are positive (repulsion between the atoms), then h is normally positive as well. However, if 51^2 are negative 
(corresponding to the less typical case of attraction between atoms belonging to a single species), then h [see Eq. (2)] 
may be either positive or negative. 

The system (1) resembles a well-known model describing the nonlinear self-phase-modulation (SPM) and cross- 
phase-modulation (XPM) interactions of light waves with different polarizations or carried by different wavelengths 
in nonlinear optics [1]. In the case of optical fibers, the evolution variable is the propagation distance z (rather than 
time), and the role of x is played by the reduced temporal variable r [1]. In optical models, however, the choice 
of the nonlinear coefficients is limited to the combinations 51 = 52 = 3/1/2 for orthogonal linear polarizations in a 
birefringent fiber and gi = 52 = h/2 for circular polarizations or different carrier wavelengths. In fact, the latter case 
is quite important in the application to two-component BECs as well, as it occurs if one assumes that the collision 
lengths for interactions between all the atoms are the same. 

Another physically interesting feature, which we include in the model to be considered below, is linear coupling 
between the two wave functions. This occurs in a mixture of two different spin states of the same isotope, which 
arises through a resonant microwave field that induces transitions between the states [5,46]. Condensates containing 
two different spin states of *^Rb have been created experimentally via sympathetic cooling [36] . In this situation, the 
normalized coupled GP equations take the form 



where the self- and cross-scattering coefficients are gi = g2 = g and h, and the linear coupling coefficient is a, which 
can always be made positive without loss of generality. 

Experimental studies of mixtures of two interconvertible condensates (with positive scattering lengths) loaded in 
an OL have not yet been reported. However, all the necessary experimental ingredients for such a work are currently 
available. Moreover, in a very recent paper [28], an experimental procedure, based on Ramsey spectroscopy and 
adjusted exactly for such a system, was elaborated. Experiments in this setting would be quite interesting, as they 
would allow the study of the direct interplay between two crucially important physical factors used as tools in the 
current experimental work — namely, the OL potential and inter-conversion between two different spin states in the 
BEC, controlled by the resonant field. Furthermore, there are now recent experimental results with linearly coupled 
BECs [23] . The use of an optical potential in the latter setting is a rather straightforward extension. 

The model combining nonlinear XPM and linear couplings, as in Eqs. (3), occurs in fiber optics as well. In that 
case, the linear coupling is generated by a twist applied to the fiber in the case of two linear polarizations, and by 
an elliptic deformation of the fiber's core in the case of circular polarizations (see, for example, [32]). However, linear 
coupling is impossible in the case of two different wavelengths. Another optical model, with only linear coupling 
between two modes, applies to dual-core nonlinear fibers, as discussed in [33] (and references therein). In the context 
of BECs, it may correspond to a special case in which the cross-scattering length is made (very close to) zero using a 
Feshbach resonance [25]. 

In this work, we aim to investigate modulated- wave states in the binary BEC described by Eqs. (3), which include 
both nonlinear and linear couplings and an OL potential. We stress that the interplay between the microwave-induced 
linear coupling in the binary model and the OL-induced periodic potential has never before been considered. As both 



dt 
dt 



(3) 
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features represent important laboratory tools, the results reported here suggest possibilities for new experiments. The 
model we study predicts new dynamical effects, such as oscillations of matter between the linearly-coupled components 
trapped in the potential wells of the OL. 

In our study, we begin by examining plane- wave solutions with V{x) = 0. When V{x) ^ 0, we apply a standing-wave 
ansatz to Eq. (3), which leads to a system of coupled parametrically forced Dufhng equations describing the spatial 
evolution of the fields. Using the method of averaging [44,42], we study periodic solutions of the latter system (called 
"modulated amplitude waves" and denoted MAWs). The stability of MAWs (and the ensuing dynamics, in the case of 
instability) is then tested by numerically simulating the underlying system of coupled GP equations. This approach, 
though simpler than the more "rigorous" computation of linear stability eigenvalues for infinitesimal perturbations, 
provides a more realistic emulation of physical experiments. Note additionally that although our stability results 
will be illustrated by a few selected examples, we have checked — by exploring different parameter regions — that these 
examples represent the MAW stability features rather generally. 

The MAW solutions are especially interesting when the system exhibits a spatial resonance. In this work, we 
consider both non-resonant solutions and solutions featuring a subharmonic resonance of the 2:1:1 form. The latter 
situation has been studied in the context of period-doubling in single-component BECs in an OL potential [30,38,39], 
but — to the best of our knowledge — spatial-resonance states in models of composite BECs have not been considered 
previously. 

An alternative (but less general) approach to the study of binary BECs with linear coupling, loaded into an OL, 
would be to seek exact elliptic-function solutions to Eqs. (3) for the case of elliptic-function potentials, V{x) = 
—Vosn^{KX,k), as has been done earlier in the two-component model without linear coupling [18]. In that work, 
stable standing-wave solutions were found under the assumption that the interaction matrix is positive definite. This 
occurs, for instance, when all the interactions are repulsive, although small negative cross-interactions are compatible 
with this condition as well. 

The rest of this paper is structured as follows: In section II, we derive plane-wave solutions and analyze their 
stability. In section III, we introduce modulated amplitude waves, and in section IV, we derive and solve averaged 
equations that describe them in both non-resonant and resonant situations. We corroborate our results and test the 
stability of the MAWs using numerical simulations. In section V, we briefly examine a more general model of a ternary 
(three-component) BEC with linear couplings. Finally, we summarize our results in section VI. 



For the linearly coupled GP equation (3), it is necessary that ki = k2 = k and /ki = /i2 = Without linear coupling, 
[as in Eqs. (1)], one may seek a broader class of solutions with independent frequencies (which correspond to chemical 
potentials in the physical context of BECs). It is important to note that the results obtained in the study of optical 
models suggest that the addition of linear coupling terms to a system of coupled NLS equations drastically alters the 
dynamical behavior [32]. In this section, we study the model without the OL, which will be included in subsequent 
sections. 

Inserting Eq. (4) into Eqs. (3) yields 



II. PLANE- WAVE SOLUTIONS 



In the absence of the external potential {V = 0), we find plane- wave solutions of the form 



tpj = Rj exp [i{kjX — ^jt)] , j = 1,2 . 



(4) 



fiRi = k^Ri + gRl + hRiRl + ai?2 , 
/xi?2 = k'^R2 + gRl + hRlR2 + aRi , 



(5) 
(6) 



which implies that 



[{g - h)RiR2 - a][Rl - Rl] = . 



One of the following two relations must then be satisfied: 



Ri = ±-^2 ; 



(7) 



R1R2 — 



a 



(8) 



g-h' 
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With Eq. (7), a nonzero solution satisfies R2 = ±-y/(fc2 + a) / {g + h). It exists when g + h > 0, provided 
k"^ + fj,^ a > 0, and when g + h < 0, if fc^ + /x =F a < 0. When g = —h, one obtains solutions of the form 
{Ri,R2) = (±-R, R) with arbitrary R and = fj,^ a. 

From Eq. (8), one finds that 




R'2 = ^±-.\ r—^ -7-^, (9) 



under the restriction that this expression must be positive. When a 7^ 0, the term inside the square root is smaller 
in magnitude than the one outside, so solutions of this type exist as long as (/x — k^) g > and the argument of the 
square root in (9) is non-negative. Hence, for repulsive and attractive BECs, respectively, the first condition implies, 
IJ,> P and fi < fc^. When h = 0, Eq. (9) takes the form 



Rl = i2g)-' (/. - fc2 ± V(/x-fc^)^-4a2) . 



For both h = and h = 2g, the condition on the argument of the square root implies that to obtain real solutions, it 
is necessary to impose the condition |/i — > 2a. 

To examine the stability of the plane waves, we substitute 

tpj {x, t) = {x, t) [1 + {x, t)], |ej 1^ <Cl , 

into Eqs. (3). This yields coupled linearized equations for ei{x,t) and e2{x,t). Assuming that Cj is periodic in x, it 
can be expanded in a Fourier series, 

00 

ej(a;,i)= ^ ejn{t) ex.^{ivnx) , 

n— — 00 

where the nth mode has wavenumber The perturbation growth rates that determine the stability of the nth mode 
are then given by 



A„ - -2iku,, ± un^-ul - .g(|i?iP + W) ± ^g^{\Ri?-\R2\^f+^h^\Ri?\R2W (10) 

where the two sign combinations ± are independent (so there are four distinct eigenvalues). Instability occurs when 
the expression under the square root in Eq. (10) has a positive real part, causing the side-band modes fc + k — Un 
of the perturbed solution to grow exponentially. In single-component condensates, this can occur only for .9 < [49]. 

Eigenvalues whose interior square root in Eq. (10) has a + sign will produce the instability before ones with a — 
sign, so we only need to check the former case. For example, if /i = 2g, the instability occurs if 



r.2+j|i?,|2 + |i?,|2_^;^;^77^^i;^i^:^j <o. (n) 



Stability conditions for the plane-wave solutions to Eqs. (3) with V = can be obtained for all the possible sign com- 
binations of g and h. We do not display them here, as they are rather cumbersome to write (although straightforward 
to compute). 

III. MODULATED AMPLITUDE WAVES 

We now generalize the above analysis to consider the two-component GP system in the presence of an OL potential. 
Toward this aim, we introduce solutions to Eqs. (3) that describe coherent structures of the form 

ijj{x,t) = Rj{x)exp[i{9ix)- ^it)] , 3 = 1,2. (12) 

Inserting Eq. (12) into Eqs. (3) and equating real and imaginary parts of the resulting equations yields 

^iRl = -R'l + Ri {d' f + gRl + V{x)Ri + hRlRi + ai?2 , 

Ati?2 = -i?2 + R2 {O' f + gR^ + V{x)R2 + hRlR2 + aRi , (13) 
O = Rie" + R[0', 

= R2e" + R'2e', (14) 
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where the prime stands for d/dx. Equations (14) imply that 

/!/ X f dx' f dx' 

with arbitrary integration constants Ci and C2, so Ri{x) = bR2{x) for some constant b unless Ci = C2 = 0. In other 
words, Ri{x) and R2{x) are different in this context only when one considers solutions with null "angular momenta." 
In the latter situation, Eqs. (13) assume the form 

R'l ~ Si , 

S[ = -ijRi + gRl + hRiRl + aR2 + V{x)Ri , 
R'2 = S2 , 

S'2 = -11R2 + gRl + hRlR2 + aRi + V{x)R2 . (16) 
When the potential V{x) is sinusoidal, Eqs. (16) are (linearly and nonlinearly) coupled cubic Mathieu equations. 



IV. AVERAGED EQUATIONS AND SPATIAL SUBHARMONIC RESONANCES 

To achieve some analytical understanding of the spatial resonances in linearly coupled BECs, we average equations 
(16) in the physically relevant case of the OL potential, 

V{x) =Vqcos{2kx). (17) 

Defining Vq = — eVb, g = eg, h = eh, and a = ea, Eqs. (16) may be written 

R'{ + ijlRi = -eVoRi cos(2Ka;) + egRl + ehRiRl + eai?2 , 

R2 + HR2 = -eVoR2 cos(2/ta;) + egRl + ehRlR2 + eaRi . (18) 
Assuming > 0, we insert the ansatz 

Rj{x) = Aj{x) cos [y/Jix) + Bj{x) sin {^/Jlx) , 

R!j{x) = —y/]lAj{x) sin [yfjlx) + ^/JlBj{x) cos {^/Jlx) (19) 

(with j = 1,2) into Eqs. (18). Differentiating the first equation of (19) and comparing it with the second yields a 
consistency condition, 

A'- cos{^x) + B'j sin(y^x) = , j = 1, 2 , 

that must be satisfied for this procedure to be valid. Inserting these equations into Eqs. (18) yields a set of coupled 

differential equations for Aj and Bj, whose right-hand sides are expanded as truncated Fourier series to isolate 
contributions from different harmonics [44,42]. The leading contribution in these equations is of 0(e), so the equations 
assume a general form 

A'j = eFA, {Ai,A2, Bx,B2,x) + 0{e^) , 

B'^ = eFB, {AuA2, ^1,^2, x) + 0{e^) . (20) 

When e = 0, Eqs. (18) decompose into two uncoupled harmonic oscillators. We have computed the exact functions 
Fa^ and Fsj in Eqs. (20) and provide them in the Appendix. 

Our objective is to isolate the parts of the functions Aj [x) and Bj [x) that vary slowly in comparison with the fast 
oscillations of cos{^fJlx) and sin(Y^a;) and to derive averaged equations governing their slow evolution. To commence 
averaging, we decompose Aj and Bj into the sum of slowly varying parts and small rapidly oscillating ones (which 
are written as power-series expansions in e): 

Aj = Aj + eWA, (Ai, A2,Bi,B2,x)+ 0{e''), 

Bj = Bj + eWB, {Ai,A2, BuB2, x) + 0{e'') . (21) 

Here, the generating functions Wa^ , Wb^ are chosen so as to eliminate all the rapidly oscillating terms in Eqs. (20) 
after the substitution of Eqs. (21). 

This procedure yields evolution equations for the averaged quantities Aj and Bj [44] , which we henceforth denote 
simply as Aj and Bj. (All other terms in the originally defined Aj and Bj are cancelled out by the choice of the 
generating functions.) As we shall see, the slow- flow equations so derived are different in resonant and non-resonant 
situations. 
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A. The Non-Resonant Case 



When ^ K [recall that k is half the wave number of the OL potential; see Eq. (17)], which is the non-resonant 
case, effective equations governing the slow evolution are 



^'1 = 4. 



e 



■^BMl + Bl) - \b, - \a,A,B^ - \b,{AI + 5Bl) 



A' 



e 



e 



-^{bMI + Bl) - §Bi - \a,A,B, - ^B,{Al + Wl) 
^A,{Al + Bl) + + \a^B,B, + ^Ai(3^2 ^ 



B' 



e 



3ff w .2 , o2^ , ^ ^A^B^B, 



-A2{Ai + Bi) 



h 



A2i3Al + B 



+ 0{e')., 
+ 0(e2). 
+ Oie'), 
+ Oie'). 



(22) 



In this case, the OL does not contribute to 0(e) terms, so the terms explicitly written in Eqs. (22) correspond to what 
one would obtain from coupled Duffing equations, as Eqs. (18) reduce to coupled Duffing oscillators in the absence 
of the OL potential [43] . These contributions yield the wavenumber-amplitude relations for decoupled condensates, 
[38,39] as well as modc-wavonumbcr relations produced by the coupling terms [44]. 



The non-resonant equations (22) give rise to three types of equilibria, at which A'^ = A'2 = B'^ 



B'2 = 0: the trivial 



(zero) equilibrium and those which we will call double modes and quadruple modes. These have, respectively, two 
and four nonzero amplitudes Aj, Bj. Single-mode and triple- mode equilibria do not exist. Different double modes 
that can be found are 7r/2 phase shifts of each other: these are "^1^2" equilibria with Ai, A2 and B\ = B2 = 0, 
and "S1B2" ones with Ai = A2 = and Bi, B2 0. 
The A-1A2 equilibria satisfy 



^^-^^-^3(5 + /^) 



(23) 



where the signs — and -|- arise, respectively, when {g + h) < 0, and {g + h) > (recall that a > 0). In the former and 
latter cases, we find that Ai = A2 and Ai = —A2, respectively. This yields the following two ^1^2 equilibria: 



(^1,^2,51,^2) =± 

(^1,^2,^1,^2) =± 



4a 



4a 



3{g + h)' V3{g + h) 



.0,0 



-4a 



-4a 



S{g + h)'V Sig + h) 



,0,0 



if g + h> 0; 
if g + h<0. 



(24) 



Similar expressions for the S1-B2 equilibria are obtained by phase-shifting the Ai^2 modes by 7r/2. 

We have examined the stability of the approximate stationary solutions corresponding to the double-mode equilibria 
obtained above with direct simulations of the coupled GP equations (3). Typically, the simulations generate solutions 
that oscillate in time (as the initial configurations are not exact stationary solutions) without developing any apparent 
instability. 

One can also find four sets of quadruple- mode equilibria in which Af = A| and Bf = B^- In the first two sets, A2 
is arbitrary: 



{A„A2,B„B2) = ±l-A2,A2,±J-Al + ^^^^^ 



{AuA2,BuB2)=± A2,^2,±4 



.2 _ 4a 
' 3{g + h) 



3{glh) 



j , if g + h>Q, 
, if g + h<0. 



(25) 
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In the second two sets, B2 is arbitrary: 



Each of the expressions (25) and (26) includes four equilibria, as there are two possible choices of the exterior signs. 
The presence of the arbitrary amplitudes in these expressions means that the quadruple-mode stationary solutions arc 
obtained as rotations of the above double- mode ones given, respectively, by Eqs. (24) and by those same equations with 
an additional 7r/2 phase shift. Accordingly, direct simulations of Eqs. (3) starting with the approximate quadruple- 
mode stationary states reveal only oscillations but no instability growth, just as with simulations initiated by the 
approximate dual-mode stationary solutions. 



B. Subharmonic Resonances 



The most fundamental spatial resonance is a subharmonic one, of type 2:1:1 [44,42,21]. In this situation, the 
parameter ^ from the initial plane-wave approximation (4) [recall that /xi = /X2 = A*] is of the form 



(27) 



where e/li is the def,uning constant [42,44,43]. [Recall that e is a small parameter; we assume /ii = 0(1).] In this 
situation, new terms occur in Eqs. (22). This leads to equations that include a contribution from the OL potential. 



A' 



B'=- 



I - Si - fsMl + Bl) \b, - \a,a,b, - \bmI + 3S|) 

- (1 + ^) ^1 + §A,{Al + Bl) + I A, + '^A.B.B, + ^A.i^Aj + Bl) 

- f f + ^) A, + |A.(^i + Bl) + + '^A.B.B, + ^A.iSAl + B^) 



+ 0{e'), 
+ 0(e2), 
+ 0{e^) , 
+ 0{e'). 



(28) 



Equations (28) have three types of equilibria when a ^ 0: the trivial one, double modes, and quadruple modes. 
When a = 0, we also find single-mode equilibria and extra double-mode ones. Triple-mode stationary solutions never 
appear. All the equilibria of Eqs. (28), except for the trivial one, correspond to spatially periodic stationary solutions 
of the underlying system (18). 

There are two kinds of A1A2 (double- mode) equilibria. The first satisfies = A2, so that the two components 
have equal amplitudes: 



{AuA2,Bi,B2) = ± 



{Ai,A2,Bi,B2) = ± 



'-4a + 2(2/zi -hFo) -4.a + 2{2iii + Vo) 



Hg + h) 



Hg + h) 



,0,0 



/4a + 2(2/xi + I/o) /4a + 2(2/xi + yo) 



3{g + h) 



Hg + h) 



,0,0 . 



(29) 



A crucial issue is the dynamical stability of these solutions, which we tested with direct simulations of the underlying 
equations (3). We found that they are stable, as exemplified in Fig. 1 for k = = 1. 
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FIG. 1. An example of evolution of the A1A2 double mode with equal amplitudes and IA2I in the case of the 2:1:1 
subharmonic resonance, constructed as per Eq. (29) for k = /x = 1, Vb = 0.1, g = h = 0.025, and a = 0.02. The left subplot 
shows the spatio-temporal evolution of IV"!!"^ by means of grayscale contour plots (|i/'2!'^ behaves similarly). The right subplot 
displays snapshots of the field {i^if for t = 99 (upper panel) and t = 105 (lower panel). In these panels, the optical-lattice 
potential is also shown by a dashed line. The results have been obtained through numerical integration of Eqs. (3) in time. 



The other A1A2 double- mode equihbrium has unequal components Ai and A2 [note that in the non-resonant case 
considered above, the double-mode equilibria, which are given by Eqs. (24) and by a n/2 phase shift thereof, always 
have equal nonzero components] : 



2 _ 2(2^1 + Vo) 



35 



> 0. 



(^1,^2,^1,-82) =± 



2fii + Vo ^ 1/(2^1+1/0)' 



16a2 



2fii + Vo 1 (2/ii + T/o)2 



16q:" 



35 



^3 



,0,0 



(30) 



where the interior -|- sign in the first component corresponds to the — sign in the second, and vice versa. The exterior 
sign ± is independent of the interior one. A necessary condition for the existence of this solution is 



2Ati + Vo 



> 



4a 



\g-h\ 



In particular, when h = 2g, which is a case of special physical relevance (as explained above), the solution becomes 

(31) 



(^l,A2,Bl,i?2) = ±( \/- 



2//1 + Vo ± \/(2/<i + Ti,)2 - 16a2 



35 



2^1 + VoT V(2mi + Vof - 16a2 



,0,0 



(32) 



provided |2/xi + Vb| > 4a. 

In fact, the existence of pairs of equilibria in which the two components have unequal amplitudes that are mirror 
images of each other is a manifestation of spontaneous symmetry breaking in the present model, which is described 
by the symmetric system of coupled equations (3). A similar phenomenon was studied in detail (in terms of soliton 
solutions) in the aforementioned model of dual-core nonlinear optical fibers, which includes only linear coupling 

between two equations [33]. 

The stability of the asymmetric stationary solutions, which correspond to equilibria with unequal components, was 
also simulated in the framework of Eqs. (3). We show the results of a typical simulation in Fig. 2. As seen in the 
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figure, these states arc subject to periodic oscillations between the two components (which is possible only in the 
presence of the linear coupling between them) . 
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FIG. 2. The same as in Fig. 1, but for the A1A2 double mode with unequal amplitudes \Ai\ and IA2I, given by Eq. (30). The 
top panel in the right subplot shows the time evolution of the numbers of atoms in the two components, iVi,2 = / \ipi,2\'^dx, 
by dashed and dash-dotted lines (the sum of the two, N = Ni + N2, is shown by the solid line). The right panel also shows 
the fields IV'il'^ and |V'2|^ (by solid and dash-dotted lines, respectively) at t = 105 and t = 210. The OL potential is shown 
by the dashed line. Oscillations of matter between the two components are clearly discernible. The parameters are g = 0.025, 
h = 0.005, a = —0.02, Vb = 0.3, and k = = 1. (Recall that the sign of a can be chosen arbitrarily.) 



The resonant equations (28) give rise to two types of B1B2 dual- mode equilibria. The first satisfies Bf = B2 and 
(^i,^2,-Bi,B2) = ± (0,0,, 



-4a + 2(2/ii - Fo) /-4a + 2(2/xi - Fq) 



3(5 + /i) 



S{9 + h) 



{Ai,A2,Bi,B2)=± \^,0, 
The second type satisfies 



/4q + 2(2/(1 - Vb) /4q + 2(2/(1 - Vo) 



3(5 + h) 



3(5 + /i) 



(33) 



Bt + Bi = 



2 _ 2(2/ii - Vo) 



35 



>0, 



{AuA2,Bi,B2) = ± 



0,0, 



' '\ 3<? 



2/ii -Vo 1 / (2fii - Vo)^ _ 16a2 



2/xi - Vo 1 {2nr-Vof 16a2 
3ff ^3V 9^ [g-hf 



(34) 



As above, the interior + sign in the first component is paired to the — sign in the second, and vice versa, whereas the 
exterior ± is independent. A necessary condition for the existence of this solution is 



2^1 - Vo 



> 



4a 



\9-h\ 



When h = 2g, the present solution becomes 



{AuA2, Bi,B2) = ±{ 0, 0, J — [2/ii -Vo± V(2mi - Vof - IGa^ 



3ff 



2^1 -VoT V(2Mi - ■^'o)^ - 16a2 



(35) 

(36) 
(37) 
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provided |2/ii — Vol > 4q:. 

Unlike the non-resonant case, the resonant S1-B2 modes are not precise phase shifts of the A1A2 modes, as the spatial 
parametric excitation resulting from the OL has only the cosine harmonic. Nevertheless, the equations describing 
these two classes of modes are similar, differing only in the sign of Vq. Direct simulations demonstrate that the 
stability of stationary solutions corresponding to the -B1-B2 equilibria is the same as in the case of the A1A2 double- 
mode equilibria considered above: the symmetric ones with = \B2\ are stable, and the asymmetric solutions with 
l-Bil ^ I-B2I are unstable. 

We have also found two sets of quadruple modes in the resonant case. The first set satisfies Bi = B2, Ai = —A2, 
and 



Al 



Vo 
2h 



35 + /1' 



a 



2h h 3g + h 



(38) 
(39) 



A necessary condition for its existence is 



> 



2h h 



3g + h 

and hence it is necessary that ^i/{3g + h) > 0. The second set of quadruple modes satisfies Bi 

Vo a , 2/11 



Ai = 



2h 



h ~^ 3g + h' 



Vo 



a 



+ i + 

A necessary existence condition for this mode to exist is 



2^1 



2/xi 



3g + h 



> 



2h 



3g + h 



-B2, Al = A2, and 
(40) 

(41) 



which also implies that Hi/{3g + /i) > 0. 

We considered quadruple modes in the presence of detuning, so /xi ^ 0. This is rather difficult to implement 
numerically, as — in view of the periodic boundary conditions in x employed in the numerical integration scheme — it 
is necessary to match both the potential and initial condition to the size of the integration domain. Nevertheless, we 
were able to perform stability simulations in this case too. We show an example of these simulations in Fig. 3. We 
observe that the quadruple mode is unstable against long-wave perturbations, even if the simulations are run with 
a = (no linear coupling). We have not observed stable quadruple states. 




10 



When a = (no linear coupling between BEC components), one can find additional double-mode equilibria and 
four single-mode ones, the latter of which take the form {Ai,0, 0, 0), (0, A2, 0, 0), (0, 0, Bi, 0), (0, 0, 0, B2), with 

A5 = ^| = Heifi±W>o, (42, 

B;=B| = -H2a^>0. (43) 

The Aj- and -Bj-modes both exist when Vo/g > 0. In this case (a = 0), matter cannot be exchanged between the 
components. In this same situation, there is also an A1B2 double-mode equilibrium [of the form {Ai,0, 0, B2)], which 
satisfies 

^'-^,-^.>'>- (") 

From Eq. (44), it follows that a necessary condition for this mode to exist is 8^1 /{3g + h) > 0. Its counterpart is an 
A2B1 equilibrium, in which the subscripts 1 and 2 are swapped in Eq. (44). 

One can extend the analysis to higher-order spatial resonances in BECs (from the lowest subharmonic resonance 
considered here) either by considering higher-order corrections to the averaged equations, or by employing a pertur- 
bative scheme based on elliptic functions, as has been done for single-component BECs in OLs [39,38]. Toward this 
aim, it may be fruitful to utilize an action-angle formulation and the elliptic-function structure of solutions to Eqs. 
(16) when Vq = 0. However, detailed consideration of higher-order resonances is beyond the scope of this work. 



V. TERNARY BECS IN OPTICAL LATTICES 



To evince the generality of the above analysis, wc briefly consider its extension to a BEC model of three hyperfine 
states coupled by two different microwave fields, which is also a physically relevant situation [19]. The corresponding 
coupled CP equations (with h = 1 and m = 1/2) are 

dip 

i-g^ = -V^Vi + S'lV'iPV'i + V{x)ipi + /ii2l'^i'2l^'^i'i + hi3\ip3\'^ipi + ai2tp2 + aisV's , 
dip 

i-g^ = -V^V2 + 5|V'2|^V'2 + V{x)lp2 + hi2\tplflp2 + /l23|'03|^V'2 + Q!l2V'l + Q!23V'3 , 

dip 

i-g^ = -V^Vs + slV'al^V's + V{x)ijj3 + /iislV'iPV's + /i23|V'2pV'3 + aisipi + a23tp2 , (45) 

where the self- and cross-scattering coefficients are g ■= gi = 92 = 93 and hjk, and the linear coupling constants are 

As in the binary case, we start with the general form (12) for stationary solutions, with 6i{x) = 62{x) = ^3(0;) = 0{x) 
and = ^2 — = Then, as done above, we set Cj = (i.e., ^ = 0) in Eqs. (15) to consider standing wave 
solutions and arrive at the following equations: 

R[ = Si , 

S[ = -uRi + gRl + hi2RiRl + hi3RiRl + OL12R2 + ai3Rz + V{x)Ri , 
R2 — S2 , 

S'2 = -11R2 + gRl + hi2RlR2 + h23R2Rl + ai2i?i + ^23^3 + V{x)R2 , 
R'3 = S3, 

S'3 = -Mi?3 + 5^3 + hl3RlR3 + h23RlR3 + ai3i?l + a23R2 + V{x)R3 , (46) 

where V{x) is the sinusoidal OL potential, as before. 

One can average Eqs. (46) with the same procedure that we applied to Eqs. (16) and thereby derive both resonant 
and non-resonant equations describing the system's slow dynamics. In particular, for the most fundamental resonant 
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case (the lowest-order, 2:1:1:1, resonance), the averaged equations are 



B' = 




4 4 



A'=- 



hi2 



/l23 



T - T - f ^^(^^ + - - - 'fA^A.B, - ■^MA.B, 



■B,iAl+3Bl)-^B,{Al + 3Bl) 



+ 0(e'), 



A'=- 



Y - T 1 ^3 - f S3(^^ + ^3^) - ^B, - ^B, - ■^A.AsB, - ■-^A.AsB, 



3ff 
8 



Q!13 o "23 

i>i 

2 2 



hi3 
4 



4 



/lis 



-■-^BsiAl + 3Bl) - ^B,{Al + 3B|) 



f + ^ Ui + |ai [Al + B?) + + ^A, + ^ + ^A,B,B, 



h 



B'=- 



- (I + ^ ) ^2 + f + B.^) + ^ + ^A, + ^AB,B, + ^A,B,Bs 



-^A,{3Ai + Bl) + ^A^i^Al + Bl) 



+ 0(e2) 



2 A2 + f^A.B.Bs + f^A^B^Bs 



.^^3(3^2 ^ ^2) ^ l^AsiiAi + Bl) 



+ 0(e 



(47) 



One can find double- mode solutions to (47) that are analogous to those of Eq. (28). For example, if \a\z \ = |a23|, 
so that the first and second components in the ternary condensate have the same strength in their linear coupling 
to the third component, there exists a double- mode equilibrium with A\ = A2 and A3 = Bi = B2 = B3 = 0. The 
values of Ai and A2 arc exactly as for binary BECs [see Eq. (29)]. except that a and h in the solution arc replaced 
by ai2 and hi2- Further, in this case, one finds Ai = —A2 for 0:13 = a23 and Ai = A2 for ais = —023. In fact, these 
modes are a straightforward extension of their two-component counterparts, as the third component is absent in the 
stationary solution. Furthermore, the stability of the symmetric double-mode equilibria, reported above, ensures the 
stability of these solutions in the ternary model. 

The situation is more interesting for asymmetric two-mode solutions, such as the ones corresponding to Eq. (30), 
which arc, simultaneously, solutions to Eqs. (45) with A3 = B3 = 0, provided Ai = — (Q2;i/«13)^2- [Note that this 
relation is used to determine 0:23 /aia, as Ai and A2 are determined from Eq. (30).] Direct simulations of the three- 
component GP equations (45) with h := h\2 = h\3 = /123 show that these asymmetric solutions are unstable, just as 
in the two-component model. The instability development, illustrated by Fig. 4. leads to an interesting dynamical 
interplay between the components. In particular, as a result of the instability, the third component is eventually 
excited, which leads to periodic oscillation of matter between all three components; i.e., in this case, we observe a 
true example of three-component dynamics. 
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FIG. 4. The same as in Fig. 2, but for the three-component (ternary) model. The bottom panel of the left plot shows the 
spatio-temporal evolution of the third component (which is absent in the unperturbed, unstable two- mode solution), and the 
thick solid line on the right shows the evolution of the number of atoms in this component (top subplot). Its spatial profile is 
shown as well (middle and bottom subplots, for t = 35.1 and t = 70.2, respectively). Periodic oscillations of matter between 
all three components are evident, cf. the oscillations in the two-compoment model, displayed in Fig. 2. The parameters are 
g = 0.025, h = 0.005, Q12 = ai3 = —0.02, Vb = 0.3, and « = The quantity 023 « 0.117 is determined from the values of 

ai2 = ai3 = —0.02 (see text). 

VI. CONCLUSIONS 

In this work, we analyzed spatial structures in coupled Gross-Pitaevskii (coupled GP) equations, which include both 
nonlinear and linear interactions, in an optical-lattice (OL) potential. The model describes a BEG consisting of a 
mixture of two different hyperfine states of one atomic species, which are linearly coupled by a resonant electromagnetic 
field. In the absence of the OL, we found plane-wave solutions and examined their stability. In the presence of the 
OL, we derived a system of averaged equations to describe a spatially modulated state which is coupled to the 
periodic potential through a subharmonic resonance. We found equilibria of the latter system and examined the 
stability of the corresponding spatially periodic solutions to the coupled GP equations using direct simulations. We 
demonstrated that symmetric dual-mode resonant states with two equal amplitudes are stable, whereas asymmetric 
ones (with unequal amplitudes) are unstable, generating solutions that oscillate periodically in time. The latter type of 
dynamical behavior is only possible in the presence of linear coupling between BEG components. Wc found four-mode 
stationary solutions as well, but they are always unstable. Finally, a three-component generalization of the model 
was introduced and briefly considered. In this case, we found that the unstable asymmetric two-mode solution, with 
one component originally empty, develops time-periodic oscillations in which the initially empty component becomes 
populated. 

ACKNOWLEDGEMENTS 

We appreciate useful discussions with Bernard Deconinck, Alex Kuzmich. Alexandru Nicolin, and Richard Rand. 
P.G.K. gratefully acknowledges support from NSF-DMS-0204585, from the Epplcy Foundation for Research and from 
an NSF-GAREER award. The work of B.A.M. was supported in a part by the grant No. 8006/03 from the Israel 
Science Foundation. 

APPENDIX 

The functions Fa^ and Fb^, which appear in Eqs. (20), can be written as a sum of harmonic contributions. To 
simpliiy the notation, we write g, h, a, and Vq simply as g, h, a, and V. 
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In the non-resonant case, = Gai / \fji>, where 
0^1(^1,^2,^3,^4,2;) = 



I - + B\) - \a,A,B, - \bMI + ^Bl) 



+ 
+ 
+ 
+ 
+ 
+ 



- - ^A,{Ai + 3Bt) - -A,B,B2 - -A^i + B^) 



|Ai(3i?2 - A\) + ^Aii?iB2 + \a^{BI - Al) 



sin(2y^a;) 

sin(4Y^a;) 



4-. 



sin(2[K; — \/Jj\x) 

h 



sin(2[K + ^/Jj\x) 



Bf + + T^BiBl 



cos{2^/Jlx) 



'-Bl + ^ 

^B,{3Al - Bf) + \a,A^B^ + ^BMl - Bl) 



V 



Bi 



cos(2kx) + 



V 



cos(2[k — y/JI\x) + 



cos(4y^x) 

cos(2[k; + y/JI\x) , 



and = Gbi/s/JI, where 
Gbi(^1, A2, ^3,^4,2;) = 



1^2 + ^A,{Al + Bl) + ^A,B,B, + ^AMl + Bl) 



+ 
+ 
+ 
+ 
+ 
+ 



1^2 + f Si(3Af + Bl) + f^A,A2B2 + \bMI + Bl) 

^-B^{iAl - Bl) + ^^1^252 + ^lB^{Al - Bl) 



sm{2^/Jix) 
sin(4y//Ix) 



sin(2[«; + v^^;) 



V A \ V 

—Bl sin(2[K-^x)+ - — 

+ |Af + ^^lA^ cos(2VMa:) 

^-AM\ - iBl) - ^A^BiB2 + '^AiiAi - Bl) 



cos(2K;a;) + 



V 



Al 



cos(2[/t — ^/Ji]x) + 



cos(4-y7i2:) 
V 



(48) 



cos(2[«; + v^a;). (49) 



Only 0(1) [i.e., constant harmonic] terms remain after averaging. 

In the resonant case, one obtains, after averaging, an extra term depending on the periodic potential V, because 
a term that was a prefactor of a non-constant harmonic in (48) and (49) has become a coefficient in front of the 
0(1) term. Other harmonic terms are also simplified due to the resonance, but they nevertheless do not contribute 
to the averaged equations because they arc still prefactors of non-constant harmonics. The extra terms with /ii arise 
from Taylor expanding in powers of e and keeping the leading-order terms. In the resonant case, Fai = Gai / k and 
Fbi = GbJk- 

In both the resonant and non-resonant cases, the expressions for Fa2 and Fb2 are obtained by switching the 
subscripts 1 < — > 2 in the equations above. 
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